Launching RStudio in Docker

First of all, use the following command to launch Rstudio in Docker.

docker run -e PASSWORD=1234 --rm -p 8787:8787 -v /Users/bfx/Documents/BCB420:/home/rstudio/docker_bcb420 risserlin/bcb420-base-image

Load the data

Use the normalized data from assignment 1.

normalized_count_data <- read.table(file=file.path("data", "GSE111958_finalized_normalized _counts.txt"),
                                    header = TRUE,sep = "\t",
                                    stringsAsFactors = FALSE,
                                    check.names=FALSE)
kable(normalized_count_data[1:5, ], type="html")

Limma

# Create data matrix
expressionMatrix <- as.matrix(normalized_count_data[,2:7])
rownames(expressionMatrix) <- normalized_count_data$GENE_NAME
colnames(expressionMatrix) <- colnames(normalized_count_data)[2:7]
minimalSet <- ExpressionSet(assayData=expressionMatrix)

# Fit the data to the model
fit <- lmFit(minimalSet, model_design)

# Apply empircal Bayes to compute differential expression for the above described model.
fit2 <- eBayes(fit,trend=TRUE)
topfit <- topTable(fit2, 
                   coef=ncol(model_design),
                   adjust.method = "BH",
                   number = nrow(expressionMatrix))

# Merge hgnc names to topfit table
output_hits <- merge(normalized_count_data[,1],
                     topfit,
                     by.y=0,by.x=1,
                     all.y=TRUE)
# Rename the colnames
colnames(output_hits)<-c("hgnc_symbol", colnames(output_hits)[2:ncol(output_hits)])

# sort by pvalue
output_hits <- output_hits[order(output_hits$P.Value),]

kable(output_hits[1:10,],type="html")

 the number of genes pass the threshold p-value < 0.05

length(which(output_hits$P.Value < 0.05))

## [1] 1476

 the number of genes pass correction

length(which(output_hits$adj.P.Val < 0.05))

## [1] 1

Account for the rep variability to improve the results

model_design_rep <- model.matrix(
  ~ samples$rep + samples$cell_type)
kable(model_design_rep,type="html")

Fit the data to the new model

fit_rep <- lmFit(minimalSet, model_design_rep)

Apply empircal Bayes to compute differential expression for the new model

# Since my data is RNA-seq data, we set the parameter trend=TRUE
fit2_rep <- eBayes(fit_rep,trend=TRUE)

topfit_rep <- topTable(fit2_rep, 
                   coef=ncol(model_design_rep),
                   adjust.method = "BH",
                   number = nrow(expressionMatrix))
# Merge hgnc names to topfit table
output_hits_rep <- merge(normalized_count_data[,1],
                         topfit_rep,
                         by.y=0,by.x=1,
                         all.y=TRUE)
# Rename the colnames
colnames(output_hits_rep)<-c("hgnc_symbol", colnames(output_hits_rep)[2:ncol(output_hits)])

# Sort by pvalue
output_hits_rep <- output_hits_rep[order(output_hits_rep$P.Value),]

kable(output_hits_rep[1:10,],type="html")

 the number of genes pass the threshold p-value < 0.05

length(which(output_hits_rep$P.Value < 0.05))

## [1] 1865

 the number of genes pass correction

length(which(output_hits_rep$adj.P.Val < 0.05))

## [1] 2

Draw MA plot, The point looks too few, but I checked dim(topfit) which is 11433, 6. Seems like the range of the data is too large.

differential_expression_results <- topfit[rownames(minimalSet), ]
differential_expression_status <- apply(differential_expression_results, 1, function(gene){
  if (gene["P.Value"] < 0.05) {
    if (gene["logFC"] < 0 ){
       return("Downregulated genes")
    }
    else{
      return("Upregulated genes")
    }
  } else {
    return("Not differentially expressed")
  }
})

limma::plotMD(fit2, column = ncol(fit2), cex = 0.5, status = differential_expression_status,
              main = "Differential gene expression")

Compare the results from two different models

simple_model_pvalues <- data.frame(hgnc_symbol = output_hits$hgnc_symbol,
                                   simple_pvalue=output_hits$P.Value)
rep_model_pvalues <-  data.frame(hgnc_symbol = output_hits_rep$hgnc_symbol,
                                 rep_pvalue = output_hits_rep$P.Value)
two_models_pvalues <- merge(simple_model_pvalues,
                            rep_model_pvalues,by.x=1,by.y=1)
two_models_pvalues$colour <- "black"
two_models_pvalues$colour[two_models_pvalues$simple_pvalue<0.05] <- "orange"
two_models_pvalues$colour[two_models_pvalues$rep_pvalue<0.05] <- "blue"
two_models_pvalues$colour[two_models_pvalues$simple_pvalue<0.05 & two_models_pvalues$rep_pvalue<0.05] <- "red"

# plot
plot(two_models_pvalues$simple_pvalue,two_models_pvalues$rep_pvalue,
     col = two_models_pvalues$colour,
     xlab = "simple model p-values",
     ylab ="Rep model p-values", 
     main="Simple vs Rep Limma")

For the interest gene

hgnc_symbol_of_interest <- "LY6E"

# There are too many overlap around the gene of interest, so I set others genes be transparent
two_models_pvalues$colour <- alpha("grey", 0.2)
two_models_pvalues$colour[two_models_pvalues$hgnc_symbol==hgnc_symbol_of_interest] <- "red"

plot(two_models_pvalues$simple_pvalue,two_models_pvalues$rep_pvalue,
     col = two_models_pvalues$colour,
     xlab = "Simple model p-values",
     ylab ="Rep model p-values",
      main="Simple vs Rep Limma (for interest gene)")

# Since there are too many overlap arround the gene of interest, use a point to mark the location of it.
points(two_models_pvalues[two_models_pvalues$hgnc_symbol==hgnc_symbol_of_interest, 2:3], pch=24,  col="red", cex=1.5)

Draw the Heatmap with the top hits using Limma (accounting for rep variability) which only show the p-value < 0.05 and columns ordered by cell type.

top_hits <- output_hits_rep$hgnc_symbol[output_hits_rep$P.Value<0.01]
heatmap_matrix_tophits <- t(
  scale(t(heatmap_matrix[which(rownames(heatmap_matrix) %in% top_hits),]))) 
heatmap_matrix_tophits<- heatmap_matrix_tophits[,
       c(grep(colnames(heatmap_matrix_tophits), pattern = "^SEU"), 
         grep(colnames(heatmap_matrix_tophits), pattern = "^SLU"))]
if(min(heatmap_matrix_tophits) == 0){
    heatmap_col = colorRamp2(c( 0, max(heatmap_matrix_tophits)), 
                             c( "white", "red"))
  } else {
    heatmap_col = colorRamp2(c(min(heatmap_matrix_tophits), 0, max(heatmap_matrix_tophits)), c("blue", "white", "red"))
  }
current_heatmap <- Heatmap(as.matrix(heatmap_matrix_tophits),
                  cluster_rows = TRUE,  show_row_dend = TRUE,
                  cluster_columns = FALSE,show_column_dend = FALSE,
                  col=heatmap_col,show_column_names = TRUE, 
                  show_row_names = FALSE,show_heatmap_legend = TRUE)
current_heatmap

 Obviously, the heatmap show that the gene expression has significant different for different cell type.

edgeR

filtered_data_matrix <- as.matrix(normalized_count_data[,2:7])
rownames(filtered_data_matrix) <- normalized_count_data$GENE_NAME

# Set up edgeR objects
d = DGEList(counts=filtered_data_matrix, group=samples$cell_type)

# Estimate Dispersion
d <- estimateDisp(d, model_design_rep)

# Fit the model
fit <- glmQLFit(d, model_design_rep)

# Calculate differential expression using the Quasi liklihood model
qlf.pos_vs_neg <- glmQLFTest(fit, coef='samples$cell_typeSLU')

kable(topTags(qlf.pos_vs_neg), type="html")

# get all results
qlf_output_hits <- topTags(qlf.pos_vs_neg,sort.by = "PValue",
                           n = nrow(normalized_count_data))

 the number of genes pass the threshold p-value < 0.05

length(which(qlf_output_hits$table$PValue < 0.05))

## [1] 498

 the number of genes pass correction

length(which(qlf_output_hits$table$FDR  < 0.05))

## [1] 4

qlf_rep_model_pvalues <- data.frame(
          hgnc_symbol = rownames(qlf_output_hits$table),
          qlf_rep_pvalue=qlf_output_hits$table$PValue)
limma_rep_model_pvalues <-  data.frame(
          hgnc_symbol = output_hits_rep$hgnc_symbol,
          limma_rep_pvalue = output_hits_rep$P.Value)
two_models_pvalues <- merge(qlf_rep_model_pvalues,
                            limma_rep_model_pvalues,
                            by.x=1,by.y=1)
two_models_pvalues$colour <- "black"
two_models_pvalues$colour[two_models_pvalues$qlf_rep_pvalue<0.05] <- "orange"
two_models_pvalues$colour[two_models_pvalues$limma_rep_pvalue<0.05] <- "blue"
two_models_pvalues$colour[two_models_pvalues$qlf_rep_pvalue<0.05 & two_models_pvalues$limma_rep_pvalue<0.05] <- "red"

# plot
plot(two_models_pvalues$qlf_rep_pvalue,
     two_models_pvalues$limma_rep_pvalue,
     col = two_models_pvalues$colour,
     xlab = "QLF rep model p-values",
     ylab ="Limma rep model p-values",
     main="QLF vs Limma ")

hgnc_symbol_of_interest <- "LY6E"

two_models_pvalues$colour <- "grey"
two_models_pvalues$colour[two_models_pvalues$hgnc_symbol==hgnc_symbol_of_interest] <- "red"
plot(two_models_pvalues$qlf_rep_pvalue,
     two_models_pvalues$limma_rep_pvalue,
     col = two_models_pvalues$colour,
     xlab = "QLF patient model p-values",
     ylab ="Limma Patient model p-values",
     main="QLF vs Limma (for interest gene)")

# Since there are too many overlaps arround the gene of interest, use a point to mark the location of it.
points(two_models_pvalues[
  two_models_pvalues$hgnc_symbol==hgnc_symbol_of_interest,2:3],
       pch=24,  col="red", cex=1.5)

top_hits <- rownames(qlf_output_hits$table)[output_hits_rep$P.Value<0.05] 
heatmap_matrix_tophits <- t(
  scale(t(heatmap_matrix[which(rownames(heatmap_matrix) %in% top_hits),]))) 
heatmap_matrix_tophits<- heatmap_matrix_tophits[,
       c(grep(colnames(heatmap_matrix_tophits),pattern = "^SEU"), 
         grep(colnames(heatmap_matrix_tophits),pattern = "^SLU"))]
if(min(heatmap_matrix_tophits) == 0){
    heatmap_col = colorRamp2(c( 0, max(heatmap_matrix_tophits)), 
                             c( "white", "red"))
  } else {
    heatmap_col = colorRamp2(c(min(heatmap_matrix_tophits), 0, max(heatmap_matrix_tophits)), c("blue", "white", "red"))
  }
current_heatmap <- Heatmap(as.matrix(heatmap_matrix_tophits),
                           cluster_rows = TRUE,
                           cluster_columns = FALSE,
                               show_row_dend = TRUE,
                               show_column_dend = FALSE,
                               col=heatmap_col,
                               show_column_names = TRUE, 
                               show_row_names = FALSE,
                               show_heatmap_legend = TRUE,
                               )

current_heatmap

# Merge gene names with the top hits
qlf_output_hits_withgn <- merge(normalized_count_data[,1],
                                qlf_output_hits, 
                                by.x=1, by.y = 0)
colnames(qlf_output_hits_withgn)<-c("hgnc_symbol", colnames(qlf_output_hits_withgn)[2:ncol(qlf_output_hits_withgn)])

qlf_output_hits_withgn[,"rank"] <- -log(qlf_output_hits_withgn$PValue,base =10) * sign(qlf_output_hits_withgn$logFC)
qlf_output_hits_withgn <- qlf_output_hits_withgn[order(qlf_output_hits_withgn$rank),]

upregulated_genes <- qlf_output_hits_withgn$hgnc_symbol[
  which(qlf_output_hits_withgn$PValue < 0.05 
             & qlf_output_hits_withgn$logFC > 0)]

downregulated_genes <- qlf_output_hits_withgn$hgnc_symbol[
  which(qlf_output_hits_withgn$PValue < 0.05 
             & qlf_output_hits_withgn$logFC < 0)]

all_differentially_expressed_genes <- qlf_output_hits_withgn$hgnc_symbol[
  which(qlf_output_hits_withgn$PValue < 0.05)]

All differentially expressed genes

The results of differentially expressed genes
The number of genesets for each data source:
GO:BP: 954
REAC: 239
WP: 16
More details for each data source:

Upregulated genes

The results of upregulated genes

The number of genesets for each data source:
GO:BP: 83
REAC: 12
WP: 4
More details for each data source:
GO:BP

REAC

WP

Downregulated genes

The results of downregulated genes

The number of genesets for each data source:
GO:BP: 918
REAC: 232
WP: 22
More details for each data source:
GO:BP

down BP

REAC

WP